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ABSTRACT 


This report presents, in two parts, a dynamic aeroelastic stability (flutter) analysis of 
a cascade of blades in supersonic axial flow. Each blade of the cascade is modeled as 
a typical section having pitching and plunging degrees of freedom. Aerodynamic 
forces are obtained from a time accurate, unsteady, two-dimensional cascade solver 
based on the Euler equations. The solver uses a time marching flux-difference 
splitting (FDS) scheme. Flutter stability is analyzed in the frequency domain. The 
unsteady force coefficients required in the analysis are obtained by harmonically 
oscillating (HO) the blades for a given flow condition, oscillation frequency, and 
interblade phase angle. The calculated time history of the forces is then Fourier 
decomposed to give the required unsteady force coefficients. An influence coefficient 
(IC) method and a pulse response (PR) method are also implemented to reduce the 
computational time for the calculation of the unsteady force coefficients for any phase 
angle and oscillation frequency. Part 1, this report, presents these analysis methods, 
and their validation by comparison with results obtained from linear theory for a 
selected flat plate cascade geometry. A typical calculation for a rotor airfoil is also 
included to show the applicability of the present solver for airfoil configurations. The 
predicted unsteady aerodynamic forces for a selected flat plate cascade geometry, and 
flow conditions correlated well with those obtained from linear theory for different 
interblade phase angles and oscillation frequencies. All the three methods of 
predicting unsteady force coefficients, namely, HO, IC and PR showed good 
correlation with each other. It was established that only a single calculation with 
four blade passages is required to calculate the aerodynamic forces for any phase 
angle for a cascade consisting of any number of blades, for any value of the oscillation 
frequency. Flutter results, including mistuning effects, for a cascade of stator airfoils 
are presented in Part 2 of the report. 
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NOMENCLATURE 


a sonic velocity 

[A j matrix of frequency domain aerodynamic coefficients 
A Roe matrix 

b airfoil semi-chord 

c airfoil chord 

ci lift coefficient 

c m moment coefficient about elastic axis 
C p pressure coefficient = 2 ( p -p^ ) / p^V^ 2 

C p steady pressure coefficient 
A C p unsteady pressure difference coefficient 
= ( P lower ~ P upper) / PooVoo 2 x amplitude 
} blade aerodynamic load vector 
F flux vector in £ direction 

G flux vector in q direction 

h plunging displacement, normal to airfoil chord 

i imaginary unit, iPT ; also index in £ direction 

i incidence angle 

[/ ] identity matrix 

I a moment of inertia about elastic axis 

Im{ } imaginary part of { } 

j index in rj direction 

kb reduced frequency, kb = cob / M cu 

[. K ] blade stiffness matrix 

Kb spring constant for plunging 

K a spring constant for pitching 

Ihh > la h > lha > la a 

frequency domain unsteady aerodynamic coefficients, Eq. (14) 
m mass of typical section 

M Mach number 

[M } blade mass matrix 

n blade index 

N number of blades in the cascade 

NB minimum number of blocks / grids required in computations 
Nk h number of reduced frequencies used in the analysis 
N(j number of interblade phase angles used in the analysis 
q vector of dependent variables 

Q transformed vector of dependent variables 

Qh hft 

Q a moment about elastic axis 

r index for interblade phase angle, Eq. (12a) 

r a radius of gyration about elastic axis in semi-chord units 

Re{ } real part of { } 

s cascade spacing or gap 

s/c gap-to-chord ratio 

S a static unbalance 
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t time 

T 2 total computational time required for two blocks calculations 
t/c thickness-to-chord ratio 

x, y Cartesian coordinates 

x a distance between elastic axis and center of mass 
in semi-chord units, see Fig. 1 

{ X } blade displacement vector 

V* reduced velocity, V* = M 00 a 00 /bco a 

a pitching displacement about elastic axis 

5 difference operator 

6 stagger angle 

p mass ratio 

p fluid density 

£ transformed coordinate 

Tj transformed coordinate 

a interblade phase angle 

t computational time, r = a^t / 2 b 

At time step 

ft) oscillation frequency 

(Oh uncoupled natural frequency for bending (plunging) 

c o a uncoupled natural frequency for torsion (pitching) 


subscripts and superscripts 


L, R left, right of interface 

T transpose of matrix 

o amplitude of harmonic motion 

co evaluated at far upstream conditions 
n evaluated at time level n 

(") dl()/dtZ 
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INTRODUCTION 


NASA Lewis Research Center has initiated an exploratory program to investigate 
the feasibility of a supersonic through flow (SSTF) fan (Refs. 1-3). The SSTF fan is 
expected to provide about a 10 percent decrease in specific fuel consumption, and 
about a 25 percent reduction in propulsion system weight, which would lead to a 22 
percent increase in aircraft range (Ref. 1). Other advantages of the SSTF fan include 
fewer fan stages required for a given pressure ratio, less inlet cowl and boundary 
layer bleed drag, better inlet pressure recovery, and better matching of bypass ratio 
variation to flight speed. For a safe design of the SSTF fan, analysis is required to 
check that dynamic aeroelastic instability or flutter will not occurr in the operating 
range. This report presents a flutter analysis for the SSTF fan blade. Only a two- 
dimensional (2D) unsteady aerodynamic / aeroelastic analysis has been developed 
and used, because of the complexity involved in modelling a three-dimensional (3D) 
rotor/ stator in supersonic through flow. 

In the 2D unsteady aerodynamic / aeroelastic analysis, the rotor / stator is modelled 
as a rectilinear two-dimensional cascade of airfoils with supersonic axial flow 
(Refs. 4-8); this is also referred to as a cascade with a supersonic leading edge locus. 
Earlier researchers have used linearized equations (see for example Ref. 4) 
applicable for unloaded cascades of flat plate airfoils in inviscid flow. The motion of 
the airfoils in each mode of the cascade is assumed to be simple harmonic with a 
constant phase angle between the adjacent blades. This assumption leads to an 
eigenvalue problem in the frequency domain and the stability of the system is 
determined from the eigenvalues (Ref 9). 

The above analyses neglect the effects of airfoil thickness, camber, and steady state 
angle of attack on supersonic cascade flutter. Recently, an approximate analysis, 
which includes these effects, was presented in Ref. 10. However, the Mach wave 
reflection pattern was assumed to be the same as that for a flat plate cascade. This 
means that the effect of the airfoil profile on the Mach wave pattern and its effect on 
stability is neglected. One way to rigorously include the effects of airfoil shape 
(thickness and camber), and steady state angle of attack in cascade flutter analysis is 
by using Computational Fluid Dynamics (CFD) methods. Recent advances in the 
field of CFD with regard to the development of more efficient algorithms and 
increases in computer speed and memory have made numerical studies to compute 
the complex flow field associated with a supersonic axial flow cascades possible. 

Steady aerodynamic characteristics of supersonic compressor airfoils using the 
Euler/Navier-Stokes equations have been reported in Ref. 11. Recently, results from 
two-dimensional CFD solvers for the analysis of oscillating cascades in supersonic 
axial flow have been reported in Refs. 12 and 13. The formulation in Ref 12 is based 
on the full-potential equation, and that in Ref 13 is based on the Euler equations. The 
discretized form of these equations are solved by the Newton-iteration method for the 
full potential equation and by an the alternating direction implicit (ADI) method for 
the Euler equations. The unsteady motion resulting from blade oscillation is included 
by generating a new grid for every time step in Ref. 12, and by using a deforming grid 
technique in Ref 13. These codes have been extended to include aeroelastic analysis 
both in the time and frequency domains in Refs. 14 and 15. 

The present study seeks to improve on the work in Ref. 13 in two respects: (a) the 
solver of Ref 13 requires explicit input of artificial viscosity factors to control the 
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numerical stability of the solution, (b) a C-grid has to be used with this version of the 
ADI solver, which inhibits the true representation of a zero-thickness flat plate 
airfoil geometry. In order to avoid these difficulties, an Euler solver based on flux 
difference splitting (FDS) has been developed in Ref. 16. This solver does not require 
explicit" input of artificial dissipation factors, and uses an H-grid which allows for a 
true representation of flat plate airfoils. This solver was modified to investigate 
aeroelastic stability of cascades in subsonic flow in Ref. 17. 

The objectives of the present study are (1) to extend the capability of the aeroelastic 
Euler cascade code of Ref. 17 to predict the unsteady aerodynamic forces of a cascade 
of airfoils in supersonic axial flow, (2) to implement influence coefficient and pulse 
response methods (Ref. 18) and verify these for supersonic axial flow cascades by 
comparison with other methods, and (3) to use the forces obtained by the above 
methods for flutter calculations of selected supersonic axial flow cascades. 

This study is restricted to a two-dimensional cascade model. A typical section 
representation of each of the blades of the cascade is assumed. Two-degrees-of- 
freedom, plunging and pitching, motion is considered for the analysis. A frequency 
domain analysis is used to predict flutter. The study is presented in two parts. In 
part I, this report, the formulation of various methods mentioned in the objectives 
above are presented and validated by comparing with linear theory results, Ref. 8. 
Typical flutter calculations, for a tuned rotor, are also presented for both a cascade of 
flat plates and actual SSTF rotor fan airfoils in supersonic axial flow. Extensive 
flutter calculations for a cascade of stator airfoils, including mistuning effects, are 
presented in part II. 


FORMULATION 

The aeroelastic model and the aerodynamic model are described in this section. 


Aeroelastic model 

The aeroelastic model for the cascade consists of a typical section with two degrees of 
freedom (bending and torsion) for each blade; see Fig. 1. Structural damping is not 
considered at present, although it could be easily included. The equations of motion 
for each blade of the cascade have the form: 

m h + S a a + K/ih = Qfi 

S^h + fa & + Ka& = Qct 

where h is the plunging (bending) displacement, a is the pitching (torsion) 
displacement, m is the airfoil mass, I a is the moment of inertia, S a is the static 
unbalance, Kh and K a are the spring constants for plunging and pitching, 
respectively; Qh and Q a are the aerodynamic loads (lift and moment, respectively) 
and the dots over the various terms indicate differentiation with respect to time. The 
above equations can be written in matrix notation as 


(la) 

(lb) 
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( 2 ) 


[*](*) + [*]{*) = (*•} 

where 

co\ 0 

0 ra OJa 



and a>fi = (Kf l /m) l/ 2 is the uncoupled natural frequency for bending; ®o =(Ka Ua) m 
is the uncoupled natural frequency for torsion; x a = S a /mb is the distance between 

the elastic axis and center of mass in semi-chord units; r a = (/« /mb 2 ) l/2 is the 
radius of gyration about the elastic axis in semi-chord units; and b is the airfoil semi- 
chord. The aerodynamic loads, Qh and Q a ,are obtained using the aerodynamic 
model described below. 


[M] = 


1 

9 

[K] = 




Aerodynamic model 

The aerodynamic model is based on the unsteady, two-dimensional, Euler equations. 
The equations in conservative differential form are transformed from a Cartesian 
reference frame (xy ) to a time dependent body-fitted curvilinear reference frame 
(£, 77 ). This transformation process and the ensuing numerical method are presented 
in detail in Ref. 19. Hence the following discussion merely highlights the 
development of the methodology and the reader is encouraged to consult Ref. 19 for 
more details. 


The transformed Euler equations can be written in vector form as 


dQ dF dG n 
dr dri 

(4) 

where 


Q = Jq 

(5a) 

F =J(StQ + &/+ Syg ) 

(5b) 

G =J(7i t q + ri x f + riyg ) 

(5c) 

and 


q = [p , pu , pv , e ] T 

(5d) 

f = [pu,pu 2 +p,puv,u(e+p)] T 

(5e) 

g = [ pv , puv , pv 2 +p ,v (e +p)] T 

(5f) 
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& > > Zy > Qt > Vx > and fly are the metrics, p is the fluid density, u and v are the 

velocities in Cartesian frame, e is the energy, p is the fluid pressure and J is the 
Jacobian of transformation. 

The approach taken in the present effort is based on the integration of these equations 
over a discrete set of contiguous cells (volumes) in computational space and is 
generally referred to as a finite volume method. This discretization results in the 
following expression where cell centers are denoted as i,j : 


dQ 5 iF 8 ijG . 
dr At, Arj 


(6a) 


With At, = Arj = 1 (by definition), this becomes 

— = -( + 8jG ) (6b) 

dr 

where 


5 i ( ) = ( ) i+ 1/2 ~ ( ) i- 1/2 and 8/ ( ) = ( )j + \a ~ ( )j- 1/2 

A consequence of the finite volume formulation is that components of the dependent 
vector Q within a particular cell represent average values over that cell. However, it 
is evident from the above representation of flux differences that a method is needed to 
allow these fluxes to be accurately represented at cell faces. As discussed in Ref. 19, 
the method used in the present effort is based on the one-dimensional Riemann 
solver of Roe (Ref. 20) at cell interfaces for each coordinate direction. The method 
uses as a basis the following approximate equation which represents a quasilinear 
form of a locally one-dimensional conservation law: 

— + A ( qt > Qr) — = 0 (7) 

dt dx 


where q is the untransformed dependent variable vector, and A {ql> Qr) is a constant 
matrix representative of local cell interface conditions which is constructed using so- 
called "Roe-averaged" variables. The determination of the eigensystem of A and 
knowing that the change in dependent variables across an interface is proportional to 
the right eigenvectors allows first order flux formulas to be constructed. This 
approach to extracting flowfield information from characteristically dictated 
directions is commonly referred to as flux difference splitting (FDS) and is applicable 
to multidimensional space if the assumption is made that all wave propagation 
occurs normal to a particular cell interface. To provide higher order spatial 
accuracy, a corrective flux is appended to the first order flux discussed above. In 
addition, in order to control dispersive errors commonly encountered with higher 
order schemes, so-called limiters are used to limit components of the interface flux 
resulting in total variation diminishing (TVD) schemes. All solutions presented 
herein were obtained using third order spatial accuracy in conjunction with the 
minmod limiter (Ref. 19). 
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An implicit scheme is used to integrate Eq. (5) in time and can be written as a 
general class of schemes (Ref. 21): 

/+ -<9AiM n 

l+y 


A Q n = 


+ f-^-[AQ n 


1+W 


U+V'/ 


( 8 ) 


A three point backward (<9 = 1, y = 1/2) method is used in the present work which 
results in second order temporal accuracy. As discussed in Ref. 19, all terms 
appearing in the above equation should be comprised of elements which result from 
one flux formula. However, experience has shown that computations performed 
using FDS theory on the right hand side of the above equation (i.e. in the residual 
term) and flux vector split (FVS) theory (Ref. 22) on the left hand side are superior to 
those obtained using FDS theory on both right and left sides. 

For the present case, 

M n = hi A + + hi A~ + hjB + + hjB~ (9) 


where 


A + = 


B + = 


/ ^\n 

dF 


\BQ I 


dG 


\ZQ 


A~ = 


B~ = 


'dF' 


\n 


\dQ 


dG 


(dQ j 


F + ,F ,G + and G are obtained using FVS theory and the residual given by 
R n = hiF n + hjG n 


( 10 ) 


is obtained by evaluating F and G using FDS theory. 

Because of the difficulty and cost of inverting the left hand side of Eq. (8), approximate 
factorization is used; detailed expressions can be found in Ref. 19. 


Grid 

The flow equations are solved on one or more passage-centered H-grids. Within a 
typical grid block (Fig. 2a), the lower computational boundary contains the upper 
surface of one blade in the cascade, while the upper computational boundary 
contains the lower surface of the adjacent blade. Periodic boundaries in the blade-to- 
blade direction extend upstream and downstream from the blade surfaces. The inlet 
boundary corresponds to the left computational boundary and the outflow 
corresponds to the right boundary. An algebraic grid generation scheme is used to 
generate a grid for a single blade passage. This grid is then stacked to form a cascade 
for multiple blades, see Fig. 2b. 


8 



The solver can simulate both pitching and plunging motions, either individually or 
in combination. The computational grid is deformed such that the airfoils follow a 
prescribed motion and the grid near the center of the passage remains fixed. This is 
done using weighting functions. For a given blade passage, the upper and lower 
boundaries containing the airfoil surfaces move according to an input function 
(which may be predetermined or calculated from the displacements obtained using a 
structural model). The grid for that passage is divided into seven sub-regions, as 
shown in Fig. 2c. The user can fix any portion of the grid near the center of the 
passage (region VII), where the grid is not affected by the airfoil displacement. The 
grid smoothly deforms for each time step in the regions between the fixed portion of 
the grid and the two blade surfaces (regions I to VI). The deforming portion of the 
grid uses a combination of weighting functions based on spatial and index 
magnitudes. For each i - index in the deforming computational plane (1 < i <ni ), a 
weighting function, wj, is calculated to control the amount of deformation in the 
j-direction: 


W JU 


s i,j 

s i, jend 



u= I V 

m -jbeg 


(x i,m ~ X i t m -1 ) + O' i, m ~ y i,m - 1 ) 


(11a) 


(lib) 


The values for jbeg and jend in Eqs. (11) depend on the region being deformed. 
Regions I through III us e jbeg = 1 and jend = jfixl. Likewise, regions IV through VI 
use jbeg = nj and jend = jfix2 (Fig. 2c). In the i-direction, similar weighting functions 
are used to keep the inlet and exit boundaries fixed: 


wiij = 

i - 1 
ile - 1 

1 

VI 

VI 

(11c) 

wi ij = 

1 

ile <i < ite 

(lid) 

wi ij = 

i -ite - 1 i 

ni -ite - 1 

ite +1 <i <ni 

(He) 


Before the grid is deformed, the coordinate system is translated to the desired center 
of pitching ( xtr,ytr ): 


Xij = Xij - xtr ij 


(Ilf) 


Yij = y ij - ytr ij 


(llg) 


The pitching and plunging motions of the blades are defined using the following 
equations: 


Xij = Xij cos(Aa) - Y ij sin(Aa) (llh) 

y t j = Yij cos(Aa) + Xij sin(Aa) + A h (Hi) 
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Here, A a and Ah are the respective pitching and plunging displacements. Finally, 
the new position of the grid (corresponding to the next time step) is calculated: 


x ij = X U + wi ij wj ij ( x ij - X ij ) + xtr ij 

y i,j = Y i,j + W H ,j w Ji,j(y i,j~Yi,j) + y tr ij 


(llj) 

(l lk) 


Boundary Conditions 

The computational grids used in the present study employ multiple blocks. A 
discussion pertaining to how flowfield conditions are imposed along the boundaries 
of the computational domain is presented next. 

Figure 2a represents a single grid block used for calculations with steady flows and 
flows with in-phase blade motions. The airfoil upper and lower surfaces are located 
along lines B-C and F-G, respectively where solid wall boundary conditions are 
employed. Lines A-E and D-H represent inflow and outflow boundaries, respectively. 
The inlet conditions are assumed to be uniform by specifying the flow density, 
velocity, flow angle, and pressure. The exit flow variables are extrapolated from the 
interior by using simple first-order model. Periodicity is imposed between lines A-B 
and E-F, and lines C-D and G-H. 

For harmonic blade motions with constant phase difference (interblade phase angle) 
between adjacent blades, N grid blocks (passages) are generally required; N is the 
number of blades in the cascade. The values of interblade phase angle (cr) that can 
occur are given, Ref. 23, as 

o r = 2Kr!N r = 0, 1,2, ••• ,N - 1 (12a) 

Then, the minimum number of grid blocks, NB, required for a given ais given as 

NB = smallest integer [ (360/cr ), 360/(360-<r ), N] (12b) 

For example, for <r = 180 deg, two grid blocks are required. This is shown in Fig. 2d. 
For this case, periodicity is enforced between lines A-B and I-J, and lines C-D and K- 
L, respectively. Also lines A-I and D-L become inflow and outflow boundaries, 
respectively. Continuity of flowfield variables is imposed between adjacent blocks (E-F 
and G-H) by simple injection (see Fig. 2d with regard to the use of interior and 
phantom cells of adjacent blocks). 

A similar procedure is followed for other interblade phase angles which generally 
require additional blade passages. However, in supersonic axial flow, the domain of 
the flow field can be reduced by inspecting the shock structures through the cascade; 
see Fig. 3. In supersonic axial flow, no disturbance can exist upstream of the leading 
edges of the blades and the wakes behind the blades cannot influence the flow in the 
blade passages. The disturbances generated inside the Mach cone of the leading edge 
do not influence the flow outside the cone, which means that only three blades need to 
be considered in the solutions. For example, consider the cascade geometry defined 
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in Fig. 3, and let the blades be denoted as shown, the blade 'ref' being the reference 
blade; let the flow conditions be defined such that the bow shock off the leading edge 
intersects the adjacent blade surfaces. The flow above blade 'ref + 1’ and below blade 
’re/'-r_will have no influence on the surface pressure on blade 'ref'. For this reason, 
the boundary conditions on the upper and lower boundaries, A-B,C-D and I-J, K-L in 
Fig. 2d, of the global grid (the grid containing two blocks) can be arbitrary. This 
means that for harmonic blade motions, only two blocks have to be considered in 
computation for any interblade phase angle. However, only the surface pressures on 
the reference blade are correct for the specified interblade phase angle. 


Flutter Analysis 

Assuming small amplitude oscillations, Eqs. (2), can be solved for flutter (stability) in 
frequency domain. The frequency domain approach, which is suitable only for linear 
problems, is described below for a tuned cascade. In a tuned cascade, all blades are 
identical and thus have identical values of structural parameters. 

The displacements for each blade of the cascade are assumed to be harmonic 
functions of time: 

[X \ = [X 0 \e i ^ t + na) » = 0, 1,2, ... ,N-1 ( 13 ) 

where N is the number of blades in the cascade, n is a blade index and co is the 
frequency of oscillation. Further, the aerodynamic loads are assumed to be linear 
functions of the displacements: 


Qh = npjy’u? [ lhh(h/b ) + lha (a) ] 

(14a) 

Q a = Ttpoob 4 a? [ lah(h/b) + l a a («) 1 

(14b) 


or 




where p = m/npo 0 b 2 is the mass ratio and [A] is the matrix of frequency domain 
unsteady aerodynamic coefficients: 

Ihh l ha 

lah la a . 

It should be noted that the elements of [A] are functions of the cascade geometry, flow 
conditions, co and <x 



The aeroelastic equations yield an eigenvalue problem (for each blade): 




( 15 ) 


u 



For a tuned cascade, the same eigenvalue problem is obtained for each blade, Ref. 9. 
Hence, it is sufficient to solve this problem for just one of the blades (but for each 

value of o ). The two eigenvalues obtained from the solution are generally complex. 
The real part of the eigenvalue determines the stability of the system; the system is 
said to flutter if the real part of either eigenvalue is positive. 

The frequency domain aerodynamic coefficients dhh , lah > etc.) are obtained from the 
unsteady Euler solver using one of the three methods described below. These 
methods are valid for small amplitude blade oscillations for which the unsteady 
flowfield is linearly dependent on the amplitude. A steady flowfield has to be obtained 
first for the given cascade geometry (stagger angle, 0and gap-to-chord ratio, sic) and 
specified inlet conditions (Mach number and incidence angle). Then the 
aerodynamic coefficients are calculated. 


1. Forced Harmonic Oscillation (HO) method 

Starting with the steady flowfield, the airfoil motions are specified as sinusoidal 
variations in time with a fixed interblade phase angle between adjacent blades: 


plunging: h= h Q sin( 2 kb + no ) 

pitching: oc= ccq sin( 2 kb M^x + no) 


(16) 

(17) 


n = 0, 1,2, ,NB-l 

where kb = cob / M^a^, is the reduced frequency based on airfoil semi-chord and 
x = aoot / 2b is a non-dimensional time. 

The Euler equations are marched in time, with the motion as specified above. After 
the initial transients have decayed, the lift and moment coefficients show a periodic 
variation with time. The variation of the lift and moment coefficients is recorded as a 
function of time. This time history is then decomposed into Fourier components to 
obtain the complex frequency-domain aerodynamic coefficients. This procedure must 
be done separately for each blade motion - plunging and pitching. It must then be 
repeated for each possible value of interblade phase angle. For a cascade with N 
blades, the values of interblade phase angle at which flutter can occur are given by 
Eq. (12a). 

As mentioned earlier, the forced harmonic oscillation method requires only two 
blocks for computation for any interblade phase angle. This is true irrespective of the 
number of blades in the cascade. However, the Euler equations have to be solved 
separately for each possible interblade phase angle, Eq. (12a). A method to obtain the 
coefficients for all possible interblade phase angles in a single calculation is 
described below. 


2. Influence Coefficient (IC) Method 

The influence coefficient method, Ref. 18, is based on the principle of linear 
superposition and it has been verified experimentally in Ref. 24. Briefly, the solution 
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to a linear problem is obtained by superposing the solutions to the individual 
elemental problems that comprise the original problem. Since the method is based on 
the principle of linear superposition, it is valid only for linear problems. 

The problem to be solved is one of a cascade of N blades in which each blade oscillates 
with a motion of the form sin(co t + no r ), where n is the blade index that varies from 0 
toiV-1, and o> is the interblade phase angle given by Eq. (12a). This problem is 
divided into N elemental problems. The k th elemental problem consists of the same 
cascade of N blades in which all blades, except one, are stationary and the k th blade 
oscillates with a motion of the form sin(<y£ ). The original problem and all the 
elemental problems have solutions that are harmonic functions of time. 

For the problem in which all blades oscillate with a motion of the form e i< ' cot+n(J r ), the 
forces (c/ and c m ) on the zeroth blade can be represented as Qoe icot ; Qq is complex 
valued to allow the force to lead or lag the motion. Denoting the force on the 0 th blade 
in the k th elemental problem as the influence coefficient Qo,* » the force Qo for a given 
o r is obtained by superposition 


JV-1 

ft>(<*) = X Qo,ke iia ' 


* = 0 


(18) 


Further, due to the spatial periodicity of the cascade, only the relative positions of the 
oscillating blade and the reference (zeroth) blade are important. That is, the forces 
generated on the 0 th blade due to the oscillation of the k th blade are equal to the forces 
on the I s * blade due to the oscillation of the k+l tfl blade, and so on. Thus, 


Qo,* = Q-k, o = Qn - *,o (19) 

where the periodicity of the cascade of N blades has been invoked again in the last 
step. Eq. (18) can be written in terms of influence coefficients as 


N 

Qo(Cr)= X QjV-A, Oe 

A = 1 


iko r 


(20a) 


Replacing the influence coefficients Qo,* by the coefficients QiV-*,0 means that all 
the required influence coefficients can be calculated simultaneously rather than 
separately. That is, instead of oscillating the k th blade, calculating the force history 
on the 0 th blade and then repeating for all values of k between 0 and JV-1, it is possible 
to oscillate the 0^ blade and determine the forces on all blades simultaneously. This 
means that the computational effort required for the calculation of all the influence 
coefficients can be reduced by a factor of N. 

For supersonic axial flow cascades, for which the aerodynamic domain of influence 
is contained within the adjacent blade passages, there is no influence on the 
reference blade from blades which are more than one blade (pitch) away. This 
indicates that it is sufficient to calculate the correct forces on just three blades, 
namely the reference blade and ones on either side. Therefore, Eq. (20) can be written 
as 
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Qref (Or) = Qref,ref + Qref+\,ref e l ° r + Qref-l,ref e lGr 


(20b) 


where the blade indices ref, ref+1 and ref-1 refer to the reference blade, the blade 
above it and the blade below it, respectively; see Fig. 4. In the present analysis, four 
blocks are used to calculate the unsteady forces on the three blades. Thus, using the 
IC method with four blocks in computations, the aerodynamic coefficients can be 
calculated for all interblade phase angles at a given oscillation frequency. This is 
true irrespective of the number of blades (iV) in the cascade. It should be noted that it 
is possible to calulate the forces on the three blades using two blocks instead of four; 
this requires some code modifications, not attempted here. 

The forced harmonic oscillation method and the influence coefficient method require 
that the Euler equations have to be solved separately for each oscillation frequency of 
interest. A method to obtain the coefficients for all frequencies in a single calculation 
is described below. 


3. Pulse Response (PR) method 

This method has evolved from the indicial approach that is widely used in many 
different fields. The indicial response is the response, lift or moment, to a step 
change in the given mode of motion. From the indicial response, the response for any 
arbitrary motion, specifically harmonic motion, can be calculated using Duhamel's 
superposition integral. 

Let the time-dependence of the blade motion be denoted as fit) and let the 
corresponding response be denoted as Fit). Let Fg it) denote the response 

corresponding to a unit step function, fit) = Sit). The response corresponding to an 
arbitrary motion f it) can be written using Duhamel’s superposition integral: 

Fit)= F git -O f'iOdC (21) 

Jo 

Using the above equation, the response to a harmonic motion, fit) = e lC0t , can be 
determined: 

Fit)= ( Fgit~Oicoe l0) Z dC (22) 

Jo 

Since only periodic response is desired, consider the above integral in the limit as 
t -> oo . Using a change of variable and extending the lower limit to -»=> , gives 

Fit) = i(oFgi(o)e icot (23) 


where Fg(co) is the Fourier transform of Fgit) given as 
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F S (t) e- i0)t dt 


(24) 


F 5 (co) = 



For an arbitrary motion fit) and the corresponding response Fit), the Fourier 
transform of Eq. (21) gives: 

F(co) = i(oF§ico) f(co) 


or 


icoFgico) = Fieri) / fieri) 


(25) 


where fieri) andFXft)) are the Fourier transforms of fit) andF(£), respectively. Eq. (25) 
indicates that the response to harmonic blade motion can be determined using any 
arbitrary motion f ( [t ) and the corresponding response Fit). 

Over the course of time, to avoid numerical difficulties, the step change in 
displacement has been replaced by a smooth step function and finally by a pulse. In 
the pulse motion, the blade returns to its original position after the duration of the 
pulse. This is in contrast to the step motion in which the blade positions are different 
before and after the step. The pulse motion thus allows the flowfield to return to its 
steady undisturbed state after the transients created by the pulse have decayed. The 
unsteady calculations therefore need to be carried out only long enough to ensure 
that the solution has reached its final state (the same as the initial state) within some 
specified tolerance. 

Out of the different pulse shapes investigated, Ref. 18 has suggested use of the 
following pulse which is used in the present calculations: 


fit) = 4 
fit) = 0 


'max 


exp 



\ 

1 

1 — t Umax i 


for 0 < £ < tmax 
for £ > t m ax 


(26) 

where t ma x is the duration of the pulse. The above choice makes fit) and fit) vanish at 
t = 0 and t = tmax in addition, higher derivatives also go to zero at t = t ma x • This 
ensures that there is a smooth transition to and from the undisturbed blade position. 

The pulse response method can be used in conjunction with the influence coefficient 
method to obtain the frequency domain aerodynamic coefficients for each possible 
interblade phase angle and for specified frequencies of oscillation, as follows. One 
blade in the cascade is given a transient motion of the form hit) = h 0 f it) or 

ait) = a Q fit). The calculations start with the steady solution and unsteady response 
to the pulse in either motion, plunging or pitching, is calculated until the transient 
flowfield reaches the steady flowfield within a specified tolerance. The motion as well 
as the response on all the blades (four in this case) are recorded and Fourier 
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transforms of these are calculated numerically for the frequency of interest. Using 
these transforms, the influence coefficients (Qk,o) are calculated from Eq. (25); it is to 

be noted that the harmonic response, icoFs(co ) , obtained from Eq. (25) for this case, is 
simply“the influence coefficient (Qk,o)- Eq. (20b) is then used to calculate the 
frequency domain unsteady aerodynamic coefficients for the interblade phase angle 
of interest. In this way, the coefficients can be determined for various values of 
reduced frequency by calculating the Fourier transforms for the frequency of interest 
using the same time histories. 

The computational implementation of the three methods, HO, IC and PR, is shown 
in Fig. 4. It shows the number of blocks used, and the associated blade motion. 


RESULTS 

The analysis methods presented in the previous section are applied to investigate the 
behavior of oscillating cascades in supersonic axial flow. First, unsteady pressure 
distributions are calculated for a flat plate cascade using the forced harmonic 
oscillations (HO) method, Eqs. (16) and (17); the results are compared with published 
results to help validate the Euler solver. Then, the unsteady lift and moment 
coefficients are compared with those obtained from the influence coefficient (IC) 
method and the pulse response (PR) method. Next, frequency domain flutter 
calculations are performed for a flat plate cascade and the results are verified by 
comparison with results from linear theory. Finally, flutter predictions are 
presented for a rotor airfoil cascade, oscillating in coupled pitching and plunging 
motion. 


All the calculations are performed for the supersonic axial flow rotor (cascade) 
geometry given in Ref. 8. The cascade has 58 blades (N= 58), a gap-to-chord ratio ( s/c ) 
of 0.31 1 and a stagger angle ( 6 ) of 28.0 deg. The inlet Mach number ( M <*,) is 2.61. All of 
the Euler solutions are obtained using a H-grid for each passage. A typical grid is 
shown in Fig. 2a for the rotor airfoil cascade and in Fig. 5 for the flat plate cascade. A 

91x41 (stream wise by pitch wise) grid is used both for flat plate airfoils and rotor 
airfoils; an algebraic grid generation scheme is used to generate this grid. The inlet 
boundary is located only 0.3 chordle^igths upstream from the leading edge because, 
in supersonic axial flow, no disturbance can travel upstream from the leading edge 
of the airfoil. There are 20 points in the streamwise direction in the region between 
the inlet boundary and the leading edge. There are 50 points along each surface of the 
airfoil, and 21 points beyond the trailing edge. The time step has been varied to check 
its effect on convergence and accuracy. The calculations have been performed on a 

Cray Y-MP computer. The solver takes about 5xl0~ 5 CPU seconds per time step per 
grid point per block. 


Code Validation 

1. Results from the forced harmonic method 
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A cascade of flat plate airfoils is considered. The incidence angle (i) is 0 deg. The 
reduced frequency based on semi-chord (kb) is 0.5 The cascade is oscillated in pitch 

about the 30% of chord location, with an amplitude of oscillation (ocq) of 0.1 deg. The 
equations are time marched for a given interblade phase angle, until a periodic 
solution is obtained. Four cycles of oscillation were found to be sufficient for the 
solutions to become periodic in time as shown in Fig. 6a. The results from the last 
cycle are Fourier decomposed to obtain the real and imaginary quantities. 

Figs. 6b and 6c show the predicted unsteady pressure distributions (first harmonic 
pressure difference coefficient, A C p , versus chord) for cr= 0 and 180 deg, respectively. 
These interblade phase angles require a minimum of one, and two grid blocks for 
numerical calculations, respectively. The results from linear theory (Ref. 4) are also 
plotted for comparison. For flat plate airfoils and small amplitudes of oscillation, the 
present calculations are expected to agree with the linear theory of Ref. 4. The results 
for the two interblade phase angles are in good agreement with those of Ref. 4. Note 
that some smearing of the pressure distribution occurs near the reflected Mach wave 
locations; however, the computed A C p 's pass through the mid-point of the high and 
low values of pressure at the discontinuity. Figure 6d shows the Mach contour lines 
predicted by the present solver. The locations at which the Mach waves strike the 

airfoil ( x/c = 0.52 on the upper surface and x/c = 0.81 on the upper sui’face) compare 
well with the locations given by linear theory. 

Results are now presented for the flat plate cascade undergoing plunging motions; 
these results are compared with results from linear theory. The amplitude of 
oscillation is, h 0 /c = 3xl0 -4 . The time response of lift coefficient is shown in Fig. 7a, 
indicating that the response has become periodic. Figures 7b and 7c show pressure 
distribution comparisons for <7=0 and 1 80 deg, respectively. As for the pitching 
motion cases, the agreement between the Euler predictions and linear theory is good. 
These comparisons validate the present Euler solver for the two interblade phase 
angles. 

It was mentioned earlier that two blocks are sufficient for calculations at any 
interblade phase angle in supersonic axial flow. To verify this, calculations are 
performed using two block, four block, and eight block configurations for ten 

interblade phase angles. It should be noted here that two blocks represent a - 0 and 
180 deg exactly; four blocks represent <7=0, 90, 180, 270 deg exactly; and eight blocks 
represent <7=0, 45, 90, 135, 180, 225, 270, 315 deg exactly. Figure 8a shows the plot of 
moment coefficient versus interblade phase angle. The cascade is pitching about 30% 
chord. It can be seen that all the three configurations give the same value indicating 
that for the supersonic axial flow solution, two blocks are sufficient to obtain 
coefficients for any interblade phase angle. Similar comparison for lift coefficient is 
shown in Fig. 8b. Results for plunging motion, not presented here, show the same 
trend. 

The variation of moment coefficient with interblade phase angle is compared with 
that obtained from the linear theory of Ref. 4, and from the ADI-Euler analysis of 
Ref. 13 in Fig. 9a. The cascade is pitching about 30% chord. Excellent agreement is 
seen with linear theory results. The discrepancy with those of Ref. 13 may be due to 
the fact that the flat plate airfoil used in the calculations of Ref. 13 was not infinitely 
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thin but had a finite thickness, which resulted in each airfoil having a blunt leading 
edge affecting the Mach wave structure. Figure 9b shows the lift coefficient versus 
interblade phase angle for the cascade oscillating in plunge. Again, excellent 
agreement is seen. As mentioned above, all the calculations are performed using 
only two blocks in the computation. 

Next, the calculated force coefficients for varying reduced frequencies are compared 
with those obtained from linear theory. Figure 10 shows this comparison for a 
pitching cascade. Figure 10a shows the lift and Fig. 10b shows the moment 

coefficients for <7=0 deg; Fig. 10c shows the lift and Fig. lOd shows the moment for 

o = 180 deg. The corresponding results for plunging motion are shown in 
Figs. 11a - lid. Excellent agreement is seen with the results from linear theory even 
for reduced frequencies as high as 2.0. It should be noted that although some 
smearing of pressure distributions occurred near shock locations (Fig. 6b-7c), all the 
force coefficients showed good agreement with the corresponding results from linear 
theory. 

This completes the validation of the present solver for supersonic axial flow 
calculations for harmonic motions. It should be noted that for the above cases, for 
unsteady runs, about 268 steps per cycle were required with four cycles of 

oscillations. This corresponds to a nondimensional time step of 8.98xl0~ 3 . For the 
amplitudes of oscillation considered, it was found that unsteady flowfleld was 
linearly related to the amplitude of oscillation, Fig. lie. Therefore, the unsteady 
aerodynamic pressures can be obtained from the influence coefficient method and 
pulse response method. The results obtained from these methods are compared with 
results from the harmonic oscillation method in the following sections. 

2. Results from the influence coefficient method 

As mentioned earlier, four blocks are used in the influence coefficient (IC) method, 
see Fig. 4b. The cascade parameters are again s/c = 0.311, 6= 28 deg, M oa = 2.61, i = 0 

deg and kb = 0.5. The reference blade in the cascade is given a forced harmonic 
motion in pitch about 30% of the chord, Eqs. (16); the remaining blades remain 
stationary. The amplitude of oscillation used in the calculations is a 0 = 0. 1 deg. Ten 
interblade phase angles are considered for calculation and comparison purposes. 
The time step used is the same as that in previous calculations, and the response on 

all the blades is calculated. Then, the required coefficients for all the values of a are 
obtained using Eq. (20b). The results from forced harmonic oscillation (HO) method 
are obtained separately for each interblade phase angle, each run taking about one- 
half hour of CPU time; the total time being about five CPU hours for the required ten 
interblade phase angles. With the influence coefficient (IC) method, it takes about 
one CPU hour, a saving factor of five. Figures 12a and 12b show the comparison of 
the lift and moment coefficients due to pitching motion at various interblade phase 
angles. The results show that the HO method and the IC method give identical 
results. It may be recalled that the HO method and the linear theory showed 
excellent agreement, as noted earlier. This verifies the superposition principle and 
the associated calculations. It implicitly confirms that the unsteady problem is linear 
for the amplitude of oscillation considered in the calculations. Figures 13a and 13b 
show the lift and moment coefficients due to plunging; good agreement is seen 
between the HO and IC methods. 
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Since the rotor cascade has 58 blades, the implementation of the HO method requires 
58 separate computations (runs) with two blocks to obtain aerodynamic coefficients 
for all interblade phase angles. With the IC method, only one computation with four 
blocks is required to obtain aerodynamic coefficients for all the 58 interblade phase 
angles. In fact, this is true for any rotor having any number of blades. However, the 
HO and IC methods have to be repeated for each oscillation frequency; this is avoided 
by using a combination of the influence coefficient (IC) and pulse response (PR) 
methods as seen below. 

3. Results from combined pulse response and influence coefficient method 

Next the pulse response (PR) method, combined with the IC method, is used to obtain 
unsteady aerodynamic coefficients at various frequencies and interblade phase 
angles. The cascade is oscillated in pitch about 30% chord. The calculations are 
performed for different pulse durations. For pulse durations of 134, 268 and 402 time 
steps, identical results were obtained for the frequencies studied. Again, only four 
blocks are used in the computations. As before, only the reference blade is given the 
pulse motion (Fig. 4c) and the response on all the blades is calculated. The 
calculations are continued until all responses have returned to the initial 
undisturbed values. The time histories are Fourier transformed and combined 
according to Eq. (25) to obtain the influence coefficients (Qk,o) which are combined 
according to Eq. (20b) to get the harmonic coefficients. Figures 14a and 14b show the 
comparison of lift and moment for pitching motion about mid-chord for a = 0 deg for 
various reduced frequencies. Excellent agreement can be seen between the 
predictions from the HO method and the PR+IC method. Figures 14c and 14d show 

the lift and moment coefficients for o= 180 deg; the same agreement is seen. The 
corresponding results for plunging motion, presented in Figs. 15a - 15d, show the 
same good agreement. Results for other interblade phase angles, not shown here, 
showed similar good agreement. This validates the implementation of the PR method 
in conjunction with the IC method. 

As mentioned above, using the IC and PR methods with four blocks, the 
aerodynamic coefficients for all 58 interblade phase angles, and all reduced 
frequencies of interest can be obtained, with minimal computational effort. The 
calculations presented above show the validity of the Euler solver, and the 
implementation of the HO, IC, and PR methods. These methods are subsequently 
used in the flutter calculations of the supersonic rotor. 

4. Relative computational costs for HO, IC and PR+IC methods 

The computational time required for the calculation of aerodynamic coefficients 
using each of the three methods can be estimated as follows. For the sake of 
conciseness, it is assumed that the total number of time steps in each of the 
calculations is the same. Let N G be the number of interblade phase angles (usually 
equal to N ), and let Nkb be the number of reduced frequencies for which calculations 
are to be performed. Let T 2 be the computational time for a single calculation using 
two grid blocks. Table I shows the total computational time for all phase angles and 
frequencies for each of the calculation methods. Compared to the HO method, the 
computational time for the IC method is reduced by a factor of N G / 2 and that for the 
PR+IC method by N a Nkb / 2. As an example, for 58 interblade phase angles and 4 
reduced frequencies, the saving factor is 29 using the IC method, and 116 using the 
PR+IC method. 
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Supersonic compressor rotor blade 


Results are now presented for a cascade of airfoils, corresponding to the rotor of a 
supersonic compressor. The cascade parameters are the same as for the flat plate 

cascade. However, the thickness-to-chord ratio of the airfoil is 0.05 and an incidence 

angle of i = 8 deg is used in the calculations. The grid used in the calculations is 
shown in Fig. 16a. 

Steady loading 

Figure 16b shows the convergence of the numerical scheme for steady solution. A 
critical examination of the numerical values showed that the forces converged to 
fourth decimal place after 2500 steps. Fig. 16c shows the distribution of steady 

pressure coefficient ( C p ) predicted by the present solver. It is compared with the 

results obtained from the ADI-Euler solver of Ref. 13. A good qualitative agreement is 
seen. Fig. 16d shows the steady pressure distribution comparison with that obtained 
from the code described in Ref. 11, which is based on a solution of the thin layer 
Navier-Stokes equations. The present solver predicts the shocks to be located about 
10% down stream of the position predicted by Ref. 11. The discrepancy is probably 
because viscous effects and quasi-3D effects are included in the analysis of Ref. 11, 
whereas the present analysis is a 2D inviscid one. 

Unsteady loading 

The unsteady pressure distribution is now calculated for the rotor cascade oscillating 
in pitch and plunge. Time history of the moment coefficient, for pitching motion 
about the 30% of chord location, is shown in Fig. 17a. The pressure difference 
coefficient (A C p ), real and imaginary, is shown in Fig. 17b for o= 0 deg and in 

Fig. 17c for cr= 180 deg. Results from linear theory are also shown for comparison. It 
should be noted that linear theory of Lane (Ref. 4) does not account for airfoil loading 
(thickness and camber). Accordingly, significant differences are seen between the 
pressure distributions predicted by the present analysis and linear theory; see 
Figs. 17b and 17c. From these figures, it can be noted that the locations of the shock 
are shifted considerably due to the airfoil camber and thickness. The shock from the 
leading edge impinges at a point 20% upstream of that predicted by linear theory and 
the reflected shock at a point 20% downstream. Figure 17d shows the Mach contour 
lines predicted by the present solver. A comparison with Fig. 6d shows the effect of 
airfoil shape on the reflected Mach wave pattern. Figure 18a shows the time response 
of lift for plunging motion; the unsteady pressure distribution is shown in Fig. 18b 
and Fig. 18c. Significant differences are seen between the present results and those 
from linear theory. 

The moment coefficient for pitching about the 30% of chord location is compared in 
Fig. 19 with that obtained from linear theory. This shows only a slight change in the 
value of the reduced frequency at which the imaginary moment coefficient becomes 
zero. This is attributed to the effect of airfoil thickness and angle of attack. It is 
surprising to note that the integrated force coefficients do not differ much from those 
predicted by linear theory even though the Mach reflection pattern and the associated 
unsteady pressure distributions are quite different. 
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Figures 20a and 20b show the results obtained from the influence coefficient and 
pulse response methods for the rotor airfoil. The results are compared with those 
obtained from the forced harmonic oscillation method. The comparisons show that 
even for airfoils with thickness and camber, the IC and PR methods give good 
correlation with HO method. 


Flutter Calculations 

The example considered is a tuned compressor cascade with 58 blades. The 
geometric and structural parameters correspond to the rotor configuration studied 
above. The structural model for each blade is a typical section (pitching and 
plunging). The pitching axis position is at 30% chord. This location was found to be 
critical for flutter in Ref. 8. The radius of gyration is r a = 0.588 and the elastic axis is 
located at the center of mass (x a - 0); the ratio of natural frequencies in bending and 
torsion is (Oh/a) a = 0.567; the mass ratio is jU = 456. 

The flutter calculations in frequency domain, (for a fixed Mach number) require an 
inner loop on interblade phase angle, with an outer loop on reduced frequency. Using 
the HO method, this can be done as follows. Starting with the steady flowfield, a value 
of reduced frequency (kb) is selected and the airfoils are forced in plunging motion 

with a specified interblade phase angle (a) until the flowfield becomes periodic in 
time. The lift and moment coefficients are calculated at each time step and stored; 
these are later decomposed into Fourier components to obtain the frequency domain 
aerodynamic coefficients. This procedure is repeated for pitching motion. In this way 
all four frequency domain aerodynamic coefficients corresponding to the specified 
values of reduced frequency and phase angle are determined. The eigenvalue 
problem given by Eq. (15) is then solved. The real part of the eigenvalue determines 

the stability; the system is stable for the selected values of kb and a if the real part of 
both eigenvalues is negative. The preceding steps are repeated for the remaining 
values of a. For a cascade with 58 blades, flutter can occur at any of the value of a 
given by Eq. (12a). If the cascade is found to be stable at all interblade phase angles, 
then a lower value of reduced frequency is selected. This is continued until the real 
part of one of the eigenvalues changes sign; this defines the flutter condition. 

With the successful implementation of the PR+IC method the unsteady aerodynamic 
coefficients for different kb and all o can be calculated in a smaller amount of CPU 
time. Only four blade passages are considered in the calculations. One of the blades 
is given a pulse, and the responses are recorded as mentioned earlier. Then using 
the PR+IC method, the coefficients are calculated for each frequency and for each 
phase angle as shown before, with considerable reduction in CPU time. The flutter 
stability calculations are verified by performing calculations for a flat plate cascade. 

Fig. 21 shows the root locus plot obtained for kb = 1.1. The plot shows the frequency 
and damping values for 58 interblade phase angles. These values of interblade phase 
angles, given by Eq. (12a), are 0, 6.2, 12.4, ..., 173.8, 180, 186.2, 192.4, ..., and 353.8 

deg; they are indicated in the figure by the arrow starting with o=0 deg.. Figure 21a 
shows the first mode, and Fig. 21b shows the second mode. The first mode is stable 
for all the values of a. The second mode is unstable for some values of c. It is 
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interesting to note that the damping values for the first mode are really one order less 
than for the second mode, even though the first mode damping values are always 
negative. The figures also show the comparison with the root locus plot obtained 
from linear theory. Excellent agreement is seen, validating the IC, PR methods, the 
flutter calculations, and the present Euler aeroelastic solver. 

Next, the calculations are repeated for the rotor airfoil cascade for kb- 1.1. The root 
locus plot is shown in Fig. 22. As before, the first mode (Fig. 22a) is more stable than 
the second mode (Fig. 22b). However, the following differences between the flat plate 
and rotor airfoil results are noted. For the first mode (Fig. 22a), the area enclosed by 
the root locus plot is more for the rotor as compared to the flat plate. The area 
enveloped is determined by the spread in the frequency and damping due to 
interblade coupling (cascade effects). Hence, Fig. 22a indicates that there is more 
cascade effects for the rotor configuration. Further, the values of damping are 

different for the rotor airfoil and flat plate for the same o. For the second mode 
(Fig. 22b), the area enveloped by the root locus plot is nearly the same for both the flat 
plate and rotor airfoils; however, the orientation of the plot is changed. Also the 

eigenvalues show an opposite trend in relation to the order of increasing a. The 
physical significance of these differences is not clear, except that they are due to the 
effect of steady loading due to thickness and camber. 

For the rotor airfoil cascade, calculations are repeated for different values of reduced 
frequency until the real part of the eigenvalue changes sign. The phase angle at 

flutter is a = 192.4 deg, the frequency ratio is oil(a a = 1.002 and the reduced frequency 
is kb = 1.09. The corresponding results for the flat plate cascade are a = 186.2 deg, 

6)/co a = 1.002 and kb = 1.153. This indicates that there is a marginal effect of airfoil 
shape on stability. 


CONCLUSIONS 

Aeroelastic stability (flutter) analysis methods using an Euler solver have been 
presented for supersonic axial flow cascades. The Euler solver is based on a flux 
difference split algorithm. Results presented have shown the validity of the unsteady 
Euler solver for harmonic blade oscillations, both in pitching and plunging motions. 
For harmonic oscillations, where the amplitude of oscillation is within the linear 
range, influence coefficient and pulse response methods have been developed. These 
methods have been successfully implemented and are shown to give same values of 
the aerodynamic coefficients as obtained from harmonic oscillation method. These 
methods considerably reduced the computational time for obtaining the aerodynamic 
coefficients required in a frequency domain based flutter analysis. It has been shown 
that a single calculation using four blade passages is enough to obtain the force 
coefficients for any interblade phase angle and frequency of interest, irrespective of 
the number of blades in the cascade. 

A two-degrees-of-freedom typical section model has been used for each blade of the 
cascade for flutter prediction. Each blade in the cascade has pitching and plunging 
motions. Representative results from a frequency domain flutter analysis have been 
obtained for a flat plate and for a rotor cascade configuration. Flutter results for the 
selected structural parameters of the case studied indicated that the airfoil shape has 
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marginal effect on stability. In part II, flutter boundaries for a cascade of stator blade 
airfoils are presented for various values of structural and geometric parameters. 
Such a parametric study will allow general conclusions to be drawn about the effect 
of airfoil shape on aeroelastic stability of supersonic flow cascades. 
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Table I: Computational costs for HO, IC and PR+IC methods 


Calculation method 

Total computational time 

HO method 


* N kb 


IC method 

1 

* N kb 

* 2 T 2 

PR+IC method 

1 

* 1 

* 2 T 2 
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Figure 2a: Cascade geometry and grid for one blade passage. 
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Figure 2b: Cascade geometry and grid for two blade passages. 
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Figure 2c: Sub-blocking used for deforming the grid. 
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Figure 2d: Periodic and injection boundaries used in calculations with multiple 
grid blocks. 
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(a) Harmonic (b) Influence (c) Pulse 

Oscillation (HO) Coefficient (IC) Response (PR) 

Fig. 4: Setup showing number of blocks used in computation with 
blade motion as indicated. 
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Fig. 5: Flat plate grid (91*41) 
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Fig. 6a: Moment coefficient vs. time in pitching motion, kb = 0.5. <7= 180 deg. 
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Fig. 6b: Unsteady pressure difference distribution, <x=0 deg, &fc = 0.5, pitching. 



x / c 


Fig. 6c: Unsteady pressure difference distribution, <r= 180 deg, kb = 0.5, pitching. 







Fig. 7b: Unsteady pressure difference distribution, a= 0 deg, kb = 0.5, plunging. 



Fig. 7c: Unsteady pressure difference distribution, <r= 180 deg, kb = 0.5, plunging. 
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interblade phase angle, a (deg) 

Fig. 8a: Comparison of moment coefficient with 2, 4 and 8 grid blocks, kb = 0.5, 

pitching. 
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Fig. 8b: Comparison of lift coefficient with 2, 4 and 8 grid blocks, kb = 0.5, pitching. 
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Fig. 9a: Variation of moment coefficient with phase angle, kb = 0.5, pitching. 
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Fig. 9b: Variation of lift coefficient with phase angle, kb- 0.5, plunging. 
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Fig. 10a: Variation of lift coefficient with reduced frequency, <7=0 deg, pitching. 



Fig. 10b: Variation of moment coefficient with reduced frequency, <r= 0 deg, 

pitching. 
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Fig. 10c: Variation of lift coefficient with reduced frequency, a = 1 80 deg, pitching. 
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Fig. lOd: Variation of moment coefficient with reduced frequency, <s- 180 deg, 

pitching. 
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Fig. 11a: Variation of lift coefficient with reduced frequency, <x=0 deg, plunging. 
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Fig. lib: Variation of moment coefficient with reduced frequency, a — 0 deg, 

plunging. 
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Fig. 11c: Variation of lift coefficient with reduced frequency, cr = 1 80 deg, plunging. 
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Fig. lid: Variation of moment coefficient with reduced frequency, a= 180 deg, 

plunging. 
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Fig. lie Variation of moment coefficient with amplitude showing linear 

relationship, kb = 1.0, c— 180 deg. 
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interblade phase angle, a (deg) 


Fig. 12a: Comparison of lift coefficient from HO and IC methods, kb — 0.5, pitching. 
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Fig. 12b: Comparison of moment coefficient from HO and IC methods, kb = 0.5, 

pitching. 
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Fig. 13a: Comparison of lift coefficient from HO and IC methods, kb = 0.5, plunging. 
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Fig. 13b: Comparison of moment coefficient from HO and IC methods, kb= 0.5, 

plunging. 
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Fig. 14a: Comparison of lift coefficient from HO and PR methods, <x = 0 deg, pitching. 



Fig. 14b: Comparison of moment coefficient from HO and PR methods, a = 0 deg, 

pitching. 
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Fig. 14c: Comparison of lift coefficient from HO and PR methods, a= 180 deg, 

pitching. 
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Fig. 14d: Comparison of moment coefficient from HO and PR methods, cr = 180 deg, 

pitching. 
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Fig. 15a: Comparison of lift coefficient from HO and PR methods, a — 0 deg, 

plunging. 
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Fig. 15b: Comparison of moment coefficient from HO and PR methods, cr = 0 deg, 

plunging. 
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Fig. 15c: Comparison of lift coefficient from HO and PR methods, a — 180 deg, 

plunging. 
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Fig. 15d: Comparison of moment coefficient from HO and PR methods, a= 180 deg, 

plunging. 
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Fig. 16 b: Conve: 
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Fig. 16c: Comparison of steady pressure distribution of the rotor airfoil, present 

solver vs Ref. 13. 
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Fig. 16d: Comparison of steady pressure distribution of the rotor airfoil, present 

solver vs Ref. 11. 
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Fig. 17a: Moment coefficient vs. time in pitching motion, kb = 1. 1, cr= 180 deg, 

a Q — 0. 1 deg. 



x / c 

Fig. 17b: Unsteady pressure difference distribution, er=0 deg, kb = 1.1, pitching. 
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Fig. 17c: Unsteady pressure difference distribution, cr= 180 deg, kb= 1.1, pitching. 



Fig. 17d: Instantaneous Mach contours for pitching rotor airfoil cascade 
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Fig. 18a: Lift coefficient vs. time in plunging motion, kb= 1.1, cr= 180 deg, 

h 0 /c= 0.0003. 
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Fig. 18b: Unsteady pressure difference distribution, er= 0 deg, kb = 1.1, pl un ging, 
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Fig. 19: Comparison of moment coefficient from HO method and linear theory, 

cr= 180 deg, pitching, rotor airfoil. 
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Fig. 20a: Comparison of moment coefficient from HO and IC methods, kb = 1.1, 

pitching, rotor airfoil. 
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Fig. 20b: Comparison of moment coefficient from HO and PR methods, < 7 = 180 deg, 

pitching, rotor airfoil. 
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Fig. 21b: Comparison of root locus plot from present analysis and linear theory, 
ki= 1.1, flat plate cascade, second mode; cr = 0. 6.2, 12.4, ..., 173.8, 180, 186.2, 192.4, 
and 353.8 deg in the direction indicated by the arrow. 
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